Main Paper
Fig2
A sigmoid function sig:
\[ p \sim \frac{1}{1+e^{-\kappa\times(x-\nu)}} \]
par(mar=c(4,4,1,1))
x=seq(0,1,.01)
inp=seq(0,1,.05)
clrs=colorRampPalette(c("purple4","yellow"))(length(inp))
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in 1:length(inp)) lines(x,sig(x,b=inp[i],a=4),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,length(inp),length.out=5)
legend("bottomright",legend=rev(sapply(inp[leg],function(d)as.expression(bquote(.(d))))),title=expression(nu),col=rev(clrs[leg]),lwd=2,cex=.8,bg="white")
tstp=25
stp=rev(seq(-3,3,length.out=tstp))
clrs=colorRampPalette(c("purple4","yellow"))(tstp)
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in rev(1:length(stp))) lines(x,sig(x,a=10^stp[i]),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,tstp,length.out=5)
legend("bottomright",legend=sapply(round(stp[leg]),function(d)as.expression(bquote(.(10^d)))),col=clrs[leg],lwd=2,title=expression(kappa),cex=.8,bg="white")
illustration of switch parameters
Fig3.tiff
A figure to illustrate different way to flatten the curve. We ran several simulationi in two different setup to show the differences between two situation: one when there is no social distancing, one in which there is.
Static simulations
The simulations have been run using:
library(parallel)
xsize=ysize=100
pop=generatePopulation(N=500,xsize=xsize,ysize=ysize,recovery=c(8,14)*25,speed=c(1,.2),behavior=rep(B,500))
cl <- makeForkCluster(18,outfile="")
neutralNSD=parSapply(cl,1:500,function(i){
print(i);
simu=slsirSimu(tstep=1500,p=c(1,1),di=2,i0=1,xsize=xsize,ysize=ysize,
pop=pop,
inf=9,
sat=1000,
p_i=1,
strategy="best",
ts=T,ap=F,visu=F
)$timeseries[,2]
})
save(file="neutralNSD.bin",neutralNSD)
pop[,"behavior"]=G
neutralSD=parSapply(cl,1:500,function(i){
print(i);
simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
pop=pop,
inf=9,
sat=1000,
p_i=1,
strategy="best",
ts=T,ap=F,visu=F
)$timeseries[,2]
})
save(file="neutralSD.bin",neutralSD)
neutralSD2=parSapply(cl,1:500,function(i){
print(i);
simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
pop=pop,
inf=9,
sat=1000,
p_i=1,
strategy="random",
ts=T,ap=F,visu=F
)$timeseries[,2]
})
save(file="neutralSD2.bin",neutralSD2)
stopCluster(cl)
plot figure
par(mar=c(2,2,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
legend("topright",legend=c("No social distancing","Social distancing"),fill=c(myred,mygreen),cex=.95)
mtext("time",1,1)
mtext("number of infected people",2,1)
load("../data/neutralSD.bin")
load("../data/neutralNSD.bin")
load("../data/neutralSD2.bin")
#we use precomputed data to compute a mean curve, those data have been generated using the scripts in "exec/graphPaper.R"
meanNSD=apply(neutralNSD,1,mean)
meanSD=apply(neutralSD,1,mean)
polygon(0:1500,meanSD,col=adjustcolor(mygreen,.5),border=NA)
polygon(0:1500,meanNSD,col=adjustcolor(myred,.5),border=NA)
lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)
abline(h=max(meanNSD),lty=2,col=myred)
abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)
text(bquote(paste("max infected (",I[max],") no SD")),x=1500,y=max(meanNSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") no SD")),at=which.max(meanNSD),adj=1,cex=.8)
text(bquote(paste("max infected (",I[max],") SD")),x=1500,y=max(meanSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") SD")),at=which.max(meanSD),adj=0,cex=.8)
abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)
Fig4.tiff
In this figure we present the overall results of the main simulations, with some some analysis of \(\delta\), the score of each simulation.
Run all simulations
To runn all simulation one need the script in exec
Those generate a lot of data, we then amde available the result of our simulation here:
you cann download
wget httpe://adsdas.ads
To get the results of the simulation then just:
name='../data/midCurveAllBadBestSLS/'
aa=lapply(unlist(lapply(list.dirs(name),list.files,pattern="allresults.bin",full.names=T)),function(f){load(f);return(allresults)})
allresults=do.call("rbind",aa)
#allresults=allresults[allresults$inf <.02 & allresults$sat >10,]
allresults=allresults[-which(allresults$max_infect <4),] #we remove the results where no outbreak happend
allresults$distances=(1-bdistance(allresults$time_max) + bdistance(allresults$max_infect))/2
For this Figure we also want to compare how simulations with different score relates to the case with no social learning. For that we re-run simulations. These simulations take time so we give the results in data/
Re-run best/worst results
library(parallel)
cl <- makeForkCluster(16,outfile="")
best=allresults[order(allresults$distances),][1:1000,]
worst=allresults[order(allresults$distances,decreasing=T),][1:1000,]
repetBest=parSapply(cl,1:500,function(i){
j=sample(nrow(best),1);print(paste(i,j));
simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
inf=best$inf[j],
sat=best$sat[j],
inf_r=best$inf_r[j],
sat_r=best$sat_r[j],
p_i=best$pind[j],
sl_rad=best$sl_rad[j],
strategy="best",
ts=T,ap=F,visu=F
)$timeseries[,2]
})
save(file="repetBest.bin",repetBest)
repetWorst=parSapply(cl,1:500,function(i){
j=sample(nrow(worst),1);print(paste(i,j));
simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
inf=worst$inf[j],
sat=worst$sat[j],
inf_r=worst$inf_r[j],
sat_r=worst$sat_r[j],
p_i=worst$pind[j],
sl_rad=worst$sl_rad[j],
strategy="best",
ts=T,ap=F,visu=F
)$timeseries[,2]
})
save(file="repetWorst.bin",repetWorst)
stopCluster(cl)
load("../data/repetBest.bin")
load("../data/repetWorst.bin")
meanWorst=apply(repetWorst,1,mean) #
meanBest=apply(repetBest,1,mean) #
generate Figure
We can now plot the graph of the paper
### left panel
par(mar=c(4,4,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
abline(v=150,lwd=2,col="gray")
par(xpd=NA)
text(150,-20,"t=150",col="1",pos=1)
par(xpd=F)
legend("topright",legend=c("No Learning (all NA)","No Learning (all A)","Worst Simulations","Best Simulations"),col=c(myred,mygreen),lty=c(1,1,2,2),lwd=2,cex=.95)
axis(1)
axis(2)
mtext("time",1,3)
mtext("number of infected people",2,3)
lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)
lines(meanBest,col=adjustcolor(mygreen,.9),lwd=2,lty=2)
lines(meanWorst,col=adjustcolor(myred,.9),lwd=2,lty=2)
### right panel
subresults=allresults[sample(nrow(allresults),20000),]
rank=rank(subresults$distances)
orderscol=rank
color_class=rev(brewer.pal(5,"PuBu"))
colorbest=color_class[1]
orderscol=rep(adjustcolor(color_class[5],.4),length(orderscol))
nvl=c(1000,2500,5000,10000)
for(i in (length(nvl):1)) orderscol[rank<=nvl[i]]=adjustcolor(color_class[i],.4)
#plot(subresults$time_max,subresults$max_infect,bg=orderscol,main=expression(I[max]*" wrt "*tau),pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
plot(subresults$time_max,subresults$max_infect,bg=orderscol,main="",pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
legend("topright",legend=c(paste("<",nvl),"all"),pt.bg=color_class,pch=21,col=adjustcolor(1,.5),pt.lwd=.1,title="Rank")
Fig 4: mean trajectories and scores for simulations
Fig5
For this figure we need the posterior distribution of our model given by the distance to the data as measured by our \(\delta\)
Get the posterior
n=1000
posterior=list(
best=allresults[order(allresults$distances),][1:n,],
worst=allresults[order(allresults$distances,decreasing=T),][1:n,],
time_max=allresults[order(allresults$time_max,decreasing=T),][1:n,],
time_max150=allresults[order(allresults$time_max150,decreasing=T),][1:n,],
time_max250=allresults[order(allresults$time_max250,decreasing=T),][1:n,],
wtime_max150=allresults[order(allresults$time_max150,decreasing=F),][1:n,],
wtime_max250=allresults[order(allresults$time_max250,decreasing=F),][1:n,],
max_infect=allresults[order(allresults$max_infect,decreasing=F),][1:n,],
max_infect150=allresults[order(allresults$max_infect150,decreasing=F),][1:n,],
max_infect250=allresults[order(allresults$max_infect250,decreasing=F),][1:n,],
wmax_infect150=allresults[order(allresults$max_infect150,decreasing=T),][1:n,],
wmax_infect250=allresults[order(allresults$max_infect250,decreasing=T),][1:n,]
)
library(hdrcde)
duocolrs=adjustcolor(c(myred,colorbest),.8)
currange=range=list(dim1=c(0,1),dim2=c(-1,3))
dimlab=list(dim1=expression(nu),dim2=expression(kappa))
### Cmpound marginal with 2d learning
expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
distribA=list(dim1=posterior$wmax_infect150$inf,dim2=posterior$wmax_infect150$sat),
distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
expnames=expnames, range=currange,cols=duocolrs,log="y",
main=""
)
expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
distribA=list(dim1=posterior$wmax_infect150$inf_r,dim2=posterior$wmax_infect150$sat_r),
distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
expnames=expnames, range=currange,cols=duocolrs,log="y",
main=""
)
expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
distribA=list(dim1=posterior$worst$inf,dim2=posterior$worst$sat),
distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
expnames=expnames, range=currange,cols=duocolrs,log="y",
main=""
)
expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
distribA=list(dim1=posterior$worst$inf_r,dim2=posterior$worst$sat_r),
distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
expnames=expnames, range=currange,cols=duocolrs,log="y",
main=""
)
Top row: distriution of parameter that gives us the best results when simulation reach 150 steps, bottom row: when the simulation reaches 1500. Left column: parameters to swtich from NA to A, right column to switch from NA to A